"""
目标：
    搜索局部密度峰值 (LDP-Searching)
输入：
    A: numpy.ndarray, shape (N, dim)，数据集
输出：
    index: ndarray of shape (N, k)，每个点的 k 近邻索引
    supk: int，自然邻居搜索的最大深度（supk）
    nb: ndarray of shape (N,), 每个点被其他点选作邻居的次数
    rho: ndarray of shape (N,), 每个点的局部密度值
    conn: ndarray of shape (N, N), 可选的自然邻域图权重矩阵
    local_core: ndarray of shape (N,), 每个点的局部核心点索引（噪声用 -10 标记）
    cores: list of int，所有局部核心点的索引
    cl: ndarray of shape (N,), 初始子簇标签
    cluster_number: int，簇的数量
"""
import numpy as np
from sklearn.neighbors import KDTree
import matplotlib
# 切换到 TkAgg 后端以兼容 PyCharm
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
from drawcluster2 import drawcluster2

def LDP_Searching(A, visualize=False):
    """
    搜索局部密度峰值。
    流程：
    1. 构建 KDTree，并找到每个点的 k=100 近邻。
    2. 循环统计不同层级时各点被当作邻居的次数，确定自然邻居 supk。
    3. 可选：构造自然邻域图，权重为 1/(1+距离)。
    4. 基于自然邻居计算局部密度 rho = nb / sum(distances)。
    5. 在 k 邻域内选取密度最大的点作为局部核心。
    6. 沿 local_core 链追溯到根，剔除噪声（标记为 -10）。
    7. 识别核心点并生成初始子簇。
    8. 可视化：画出多种代表点与簇结构图。

    参数：
        A: numpy.ndarray, shape (N, dim)
        build_conn: bool, 是否构造自然邻域图
        visualize: bool, 是否绘制示意图

    返回：
        如果 build_conn=True，则返回 index, supk, nb, rho, conn, local_core, cores, cl, cluster_number
        否则，返回 index, supk, nb, rho, local_core, cores, cl, cluster_number
    """
    # 转为 numpy 数组
    A = np.asarray(A)
    N, dim = A.shape

    # 1. 构建 KDTree 并查询 k=100 近邻
    tree = KDTree(A, leaf_size=1, metric='euclidean')  # leaf_size指定树中叶子节点的最大点数，选用欧式距离
    k = 100 # 经验值
    index = tree.query(A, k=k, return_distance=False) # 返回排好序的索引和距离

    # 2. 统计自然邻居次数
    nb = np.zeros(N, dtype=int) # 自然邻居个数
    r = 1
    count = 0   # 自然最近邻数为零的数据量连续相同的次数
    count1 = 0  # 前一次自然最近邻数为零的数据量
    count2 = 0  # 此次自然最近邻数为零的数据量
    while True:
        for i in range(N):
            nb[index[i, r]] += 1
            # RNN[k, nb[k]] = i 假设 RNN 是一个二维 NumPy 数组，k 和 nb 是索引
        r += 1
        # count2 = 0
        count2 = np.sum(nb == 0)
        # 计算nb(i)=0的点的数量连续不变化的次数
        count = count + 1 if count2 == count1 else 1    # 如果连续两次自然最近邻数为零的数据量相同，则计数加一，否则重置为 1
        # 邻居搜索终止条件
        if count2 == 0 or (r > 2 and count >= 2):
            break
        count1 = count2

    # 计算自然最近邻的各种特征量
    supk = r - 1  # 最终K值，也是自然最近邻居的平均数
    min_nb = np.min(nb)  # 自然邻居的最小数目
    max_nb = np.max(nb)  # 自然邻居的最大数目
    # NN = index[:, 1:supk + 1]  # 各数据点的K近邻数据点集
    # ratio_nb = nb / (N * supk)  # 各数据点的自然最近邻居数目所占比例
    print("supk:" + str(supk))
    print("min_nb:" + str(min_nb))
    print("max_nb:" + str(max_nb))

    # 初始化密度数组
    rho = np.zeros(N)
    Non = max_nb

    # 3. 构造自然邻域图
    conn = None
    if visualize:
        conn = np.zeros((N, N), dtype=float)
        for i in range(N):
            for j in range(1, supk + 1):
                nbr = index[i, j]
                dist = np.linalg.norm(A[i] - A[nbr])
                weight = 1.0 / (1.0 + dist)
                conn[i, nbr] = weight
                conn[nbr, i] = weight

    # 4. 计算局部密度 rho
    max_nb = nb.max()
    rho = np.zeros(N, dtype=float)
    for i in range(N):
        dist_sum = 0.0
        for j in range(max_nb + 1):
            dist_sum += np.linalg.norm(A[i] - A[index[i, j]])   # 计算欧几里得距离
        rho[i] = nb[i] / dist_sum if dist_sum > 0 else 0.0      # 计算密度

    # 5. 选取局部核心点（找到每个点的代表点：搜索每个点的k邻域，将其中密度最大的点作为其代表点）
    local_core = np.zeros(N, dtype=int)     # 存放n个点的局部核点
    for i in range(N):
        rep = i
        max_rho = rho[i]
        for j in range(1, supk + 1):
            neighbor = index[i, j]
            if rho[neighbor] > max_rho:
                max_rho = rho[neighbor]  # 当前密度值小于邻居的密度值，将密度值更大的邻居作为代表点
                rep = neighbor
        local_core[i] = rep     # 最终的代表点

    # 6. TODO 追溯代表链（并剔除噪声）
    # 数据集名字                 簇数      迭代次数      筛选排除的round_id
    # aggregation.txt           7         699          686
    # d6.txt                    4         1307         1283
    # E6.txt                    7         7762         7622
    # t4.txt                    6         7149         7020
    # t7.txt                    9         7139         7010
    # t8.dat                    8         7204         7074
    visited = np.zeros(N, dtype=int)
    round_id = 0
    for i in range(N):
        if visited[i] or local_core[i] == -10:
            continue    # 如果已经访问过或者是噪声点，则跳过
        round_id += 1
        parent = i
        while local_core[parent] != parent and local_core[parent] != -10:
            visited[parent] = round_id
            parent = local_core[parent]

        # TODO 剔除噪声点的 round 阈值
        # max_depth = N - 0.018 * N
        max_depth = 7190
        # new：
        if round_id > max_depth:     # 假设 round 大于 ?? 的点视为噪声点
            local_core[visited == round_id] = -10 # 将噪声点的代表点设置为 -10
        else:                   # 如果 round 较小，正常进行代表点分配
            local_core[visited == round_id] = parent # 将当前轮次的元素赋值为 parent

        # old：
        # local_core[visited == round_id] = parent # 将当前轮次的元素赋值为 parent
    print("round_id:" + str(round_id))


    # 7. 核心点识别与子簇分配
    cores = []
    cl = np.zeros(N, dtype=int)
    cluster_number = 0
    for i in range(N):
        if local_core[i] == i:
            cluster_number += 1
            cores.append(i)
            cl[i] = cluster_number
    for i in range(N):
        if local_core[i] >= 0:
            cl[i] = cl[local_core[i]]

    # 8. 可视化多图
    if visualize:
        plt.rcParams['font.sans-serif'] = ['SimHei']  # Ensure Chinese font support
        plt.rcParams['axes.unicode_minus'] = False  # Fix for negative signs
        # 图1：所有点与代表点连线
        plt.figure()
        plt.scatter(A[:, 0], A[:, 1], c='k', s=10)
        for i in range(N):
            p = local_core[i]
            if p != -10:
                plt.plot([A[i, 0], A[p, 0]], [A[i, 1], A[p, 1]], linestyle=':', color='gray', linewidth=1)
        rep_pts = [i for i in range(N) if local_core[i] == i]
        plt.scatter(A[rep_pts, 0], A[rep_pts, 1], marker='s', s=30, facecolors='none', edgecolors='r')
        plt.title('代表点连线图')
        plt.show()

        # 图2：跳过噪声点的代表点图
        plt.figure()
        valid = local_core != -10
        plt.scatter(A[valid, 0], A[valid, 1], c='k', s=10)
        for i in range(N):
            p = local_core[i]
            if p != -10:
                plt.plot([A[i, 0], A[p, 0]], [A[i, 1], A[p, 1]], linestyle=':', color='gray', linewidth=1)
        plt.scatter(A[rep_pts, 0], A[rep_pts, 1], marker='s', s=30, facecolors='none', edgecolors='r')
        plt.title('代表点（噪声已排除）')
        plt.show()

        # 图3：核心点图以及相应的初始聚类结果
        plt.figure()
        plt.scatter(A[:, 0], A[:, 1], c='lightgray', s=10)
        for i in range(N):
            p = local_core[i]
            if p >= 0:
                plt.plot([A[i, 0], A[p, 0]], [A[i, 1], A[p, 1]], linestyle='-', color='gray', linewidth=0.5)
        drawcluster2(A, cl, cluster_number + 1)
        plt.scatter(A[cores, 0], A[cores, 1], marker='s', s=50, facecolors='none', edgecolors='r')
        plt.title('搜索 supk 近邻的结果')
        plt.show()

        # 图4：自然邻域图
        N, dim = A.shape
        assert dim == 2, "只支持 2D 数据绘图"
        plt.figure()
        plt.scatter(A[:, 0], A[:, 1], c='k', s=10, label='Data Points')
        # 遍历 conn 矩阵，绘制边
        for i in range(N):
            for j in range(i + 1, N):  # 避免重复绘制对称边
                if conn[i, j] > 0:  # 仅绘制有权重的边
                    plt.plot([A[i, 0], A[j, 0]], [A[i, 1], A[j, 1]],
                             linestyle='-', color='gray', linewidth=conn[i, j] * 2)
        plt.title("Natural Neighborhood Graph")
        plt.xlabel("X")
        plt.ylabel("Y")
        plt.legend()
        plt.show()

    # 返回结果
    return index, supk, nb, rho, local_core, cores, cl, cluster_number
